#ifndef CUDA_SPH
#define CUDA_SPH

#include <glm/glm.hpp>
#include <vector>


// Interface to the outside world
struct CudaParticle {
	CudaParticle() {
		pos = glm::vec3(0,0,0);
		vel = glm::vec3(0,0,0);
		currTime = 0;
		otLev = 0;
		cntr = 0;
		lastLev = 0;
	}
	glm::vec3 pos;
	glm::vec3 vel;

	// Testing adaptive timestep
	float currTime;

	// For the octree
	// This is the level of the octree the particle is in.
	int otLev;

	int cntr;
	int lastLev;
};

class CudaSPH {
public:

	CudaSPH( glm::vec3 & gridMin, glm::ivec3 numCells, glm::vec3 cellWidth );

	void SetSPHParams( 
		float ts, float k, float visc,
		float gamma, float supportRad);

	int GetNumParticles() {
		return this->m_part_num;
	}

	// Allocates particles.
	void SetParticles(const std::vector<CudaParticle> & inParts)
	{
		m_part_num = inParts.size();

		// Allocate particles
		AllocateParticlesDevice(inParts);

		m_firstRun = true;
	}

	void SetMass( float inMass) {
		m_mass = inMass;
	}

	void AddParticles(const std::vector<CudaParticle> & inParts ){
		std::vector<CudaParticle> newParts;
		// Gets the current particles
		GetParticles(newParts);
		int numOldParts = newParts.size();
		
		newParts.resize(numOldParts + inParts.size());

		std::copy(inParts.begin(), inParts.end(), newParts.begin() + numOldParts );
		
		SetParticles(newParts);

		// We don't want to recompute the density, so we set m_firstRun to false;
		if (numOldParts != 0){
			m_firstRun = false;
		}
	}

	// Returns the particles from the device.
	void GetParticles(std::vector<CudaParticle> & outParts);

	// Returns the grid ndx 
	void GetGridNdx( std::vector<int> & outGrid);

	// Returns the grid cells
	void GetGridEnd( std::vector<int> & outGrid);

	// Returns the forces
	void GetForces(std::vector<float> & outForce );

	// Returns pressure
	void GetPress(std::vector<float> & outPress);

	// Returns density
	void GetDens(std::vector<float> & outDens);

	// Returns the octree EoF
	void GetOctreeEoF(std::vector<unsigned int> & outOctreeEoF);

	unsigned int GetOctreeNumLevs();

	void ClearOctree();

	// Takes an update step.
	void Update( );

private:
// Internal functions.
	void DeallocateParticles();

	// Gets the average of a vector
	float GetAverageSerial( const std::vector<float> & inVect );

	// Allocates the two grid vectors on the GPU.
	void AllocateGridDevice( glm::vec3 & gridMin, glm::ivec3 & numCells, glm::vec3 & cellWidth);

	// Allocates the particles on the GPU.
	void AllocateParticlesDevice(const std::vector<CudaParticle> & inParts );

	// Allocates memory for the octree
	void AllocateOctree();

// Variables
	
	//
	//// Device arrays.
	CudaParticle * mD_parts;
	int * mD_parts_gNdx;
	float * mD_parts_press;
	float * mD_parts_dens;
	float * mD_parts_forces;

	// Octree
	unsigned int * mD_octree_EoF;
	unsigned int * mD_octree_sndx;
	unsigned int * mD_octree_eend;
	unsigned int m_totOctreeCells;

	// Particle indices for the grid and the counts.
	int * mD_grid_ndx;
	int * mD_grid_end;

	//
	//// Global variables
	float m_mass;
	int m_part_num;
	bool m_firstRun; //Is this the first run or not.  we calculate rest density on the first run.

};

#endif